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SUMMARY 


Galerkin  and  least-squares  finite  element  formulations  in 
terms  of  the  primitive  variables  have  been  applied  to  the 
equations  governing  compressible,  inviscid  flow.  A novel 
finite  element  representation  for  the  groups  of  variables, 
rather  than  the  single  variables>  occurring  linearly  in  the 
conservation  form  of  the  governing  equations  has  led  to  a 
relatively  sparse  stiffness  matrix.  The  Galerkin  formulation 
was  used  in  conjunction  with  Newton's  method  but  solutions 
for  the  flow  about  circular  cylinders  could  only  be  obtained 
with  freestream  Mach  numbers  less  than"  0.32.  The  least- 
squares  formulation  was  applied  in  conjunction  with  an 
iterative  scheme  of  the  successive  over- relaxation  type. 
Solutions  have  been  obtained  for  the  flow  about  circular  and 
elliptic  cylinders,  a 6%  circular-arc  aerofoil  and  a NACA-0012 
aerofoil  at  zero  angle  of  attack,  with  the  free-stream  Mach 
number  sufficiently  large  that  locally  sonic  conditions  have 
occurred.  The  solutions  are  in  good  agreement  with  both 
experimental  results  and  other  computational 
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Galerkin  and  least-squares  finite  element  formulations  in  terms  of  the 
primitive  variables  have  been  applied  to  the  equations  governing 
compressible,  inviscid  flow.  A novel  finite  element  representation  for 
the  groups  of  variables,  rather  than  the  single  variables,  occurring 
linearly  in  the  conservation  form  of  the  governing  equations  has  led  to  a 
relatively  sparse  stiffness  matrix.  The  Galerkin  formulation  was  used  in 
conjunction  with  Newton's  method  but  solutions  for  the  flow  about  circular 
cylinders  could  only  be  obtained  with  freestream  Mach  numbers  less  than 
0.32.  The  least-squares  formulation  was  applied  in  conjunction  with  an 
iterative  scheme  of  the  successive  over- relaxation  type.  Solutions  have 
been  obtained  for  the  flow  about  circular  and  el  liptic  cylinders , a 6% 
circular-arc  aerofoil  and  a NACA-0012  aerofoil  at  zero  angle  of  attack, 
with  the  free-stream  Mach  number  sufficiently  large  that  locally  sonic 
conditions  have  occurred.  The  solutions  are  in  good  agreement  with  both 
experimental  results  and  other  computational  solutions. 
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1 . INTRODUCTION 

The  purpose  of  this  report  is  to  examine  a number  of  different  finite  element 
formulations  suitable  for  external,  inviscid,  subcritical  flow.  As  such  it  is 
a stepping-stone  towards  the  treatment  of  transonic,  external,  inviscid  flow. 

Most  previous  finite  element  applications  to  subcritical  flow  e.g. (ref. 1,2,3) 
have  been  based  on  either  a velocity  potential  or  stream  function  formulation. 
However  no  formulation  based  on  the  full  equations  of  motion  has  been  successful 
in  obtaining  solutions  that  contained  significant  regions  of  embedded  supersonic 
flow.  For  this  and  other  reasons  set  out  in  Section  2,  the  present  formulation 
makes  direct  use  of  the  primitive  variables,  i.e.  velocity,  density  and 
pressure.  It  is  believed  that  this  is  the  first  time  a primitive  variable 
finite  element  formulation  has  been  applied  to  the  full  equations  of  motion 
governing  compressible,  inviscid  flow. 

In  order  to  ensure  an  efficient  finite  element  algorithm  it  is  important  to 
establish  consistent  analytic  representations  for  each  of  the  dependent  variables; 
this  problem  is  considered  in  Section  2.  For  problems  that  do  not  possess  a 
variational  formulation  the  two  most  effective  finite  element  methods  are  based 
on  the  Galerkin  and  least-squares  formulations.  Both  these  formulations  have 
been  applied  to  the  current  problem  and  reduction  of  the  governing  partial 
differential  equations  to  algebraic  equations  is  described  in  Section  2. 

The  governing  equations,  in  both  the  differential  and  algebraic  form,  are 
highly  nonlinear  and  require  iterative  solution  techniques;  these  are  described 
in  Section  3.  A modified  Newton's  method  has  been  used  in  conjunction  with  the 
Galerkin  formulation.  Because  Newton's  method  will  only  converge  for  a starting 
solution  relatively  close  to  the  final  solution  it  has  been  necessary  to  use  an 
ancilliary  technique  to  get  close  to  the  final  solution.  This  has  been  achieved 
by  utilising  the  unsteady  version  of  the  governing  equations.  After  application 
of  the  finite  element  formulation,  these  equations  have  been  treated  as  ordinary 
differential  equations  in  time  and  have  been  integrated  until  the  converged 
solution  is  approached.  Because  the  least-squares  formulation  leads  to  a 
positive-definite  stiffness  matrix  it  has  been  possible  to  apply  a successive 
over-relaxation  (SOR)  method  of  solution  after  locally  linearising  the  algebraic 
equations . 

Solutions  for  the  flow  about  circular  and  elliptic  cylinders  using  both  the 
Galerkin  and  least-squares  formulations  are  given  in  Section  4.  Solutions  for 
the  flow  about  two  representative  aerofoils  at  various  free-stream  Mach  numbers 
are  presented  in  Section  5.  Comparisons  are  made  with  experimental  results  and 
with  other  computational  results.  A comparison  of  the  Galerkin  and  least-squares 
formulations  as  applied  to  the  current  problem  is  made  in  Section  6. 


2.  FORMULATION  OF  THE  PROBLEM 

In  Section  2.1  the  relative  merits  of  the  velocity  potential  and  primitive 
variable  formulations  applied  to  subcritical,  inviscid  flow  are  presented. 

The  equations  of  motion  appropriate  to  both  steady  and  unsteady  flow  and  the 
corresponding  boundary  conditions  are  indicated  in  Section  2.2.  Different  orders 
of  analytic  representation  for  different  variables  are  possible  and  these  are 
discussed.  Both  the  Galerkin  and  least-squares  finite  element  formulations  have 
been  used  to  reduce  the  governing  partial  differential  equations  to  algebraic 
equations;  these  are  described  in  Sections  2.3  and  2.4  respectively. 

2.1  Velocity  potential  vs.  primitive  variables 

A velocity  potential  formulation  of  subcritical,  inviscid  flow  might 
express  the  governing  equations  in  terms  of  the  velocity  potential,  \p  and  the 
local  sound  speed,  a.  In  two  dimensions  the  governing  equations  are 
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^ • V = " 

and 

a^  + = a^  . (2) 

2 '■X  y^  o ^ 

In  equations  (1)  and  (2)  y is  the  specific  heat  ratio  and  a^  is  the  stagnation 

value  of  the  sound  speed.  The  use  of  equations  (1)  and  (2)  is  attractive 
because  only  one  unknown,  <P , is  required  at  each  node.  Equation  (2)  is 
treated  as  a local  algebraic  relationship  which  adjusts  the  value  of  'a'  at 
each  step  of  the  iteration. 

Equation  (1)  can  be  written  as  a Poisson-like  equation 


XX 


+ VJ, 


yy 


= f. 


(3) 


in  which  f incorporates  all  the  terras  associated  with  the  compressible  nature 
of  the  flow.  Application  of  a variational  finite  element  formulation  of 
equation  (3j  permits  an  iterative  solution  in  which  f is  recalculated  after 
each  solution  for  <fi.  The  iterative  procedure  fails  at  a local  Mach  number 
of  unity(ref.2,3) . This  is  probably  due  to  the  use  of  equation  (2)  to 
adjust  the  value  of  'a'.  It  is  interesting  that  a linearisation  of 
equation  (1)  that  leads  to  the  transonic  small  perturbation  equation  and 
avoids  the  use  of  equation  (2)  can  be  iterated  to  locally  supersonic  Mach 
numbers (ref. 4) . 

Another  disadvantage  of  the  use  of  equations  (1)  and  (2)  is  that  the 
converged  solution,  must  be  differentiated  numerically  before  a useful 
quantity,  the  pressure  at  the  body  surface,  can  be  obtained.  Also  since 
equation  (1)  is  highly  nonlinear,  application  of  a finite  element  formulation 
results  in  a large  number  of  cross-terms  that  must  be  manipulated  at  each 
step  of  the  iterative  process.  This  results  in  a considerable  increase  in 
computation  time. 

To  avoid  some  of  the  disadvantages  noted  above  the  present  formulation 
makes  use  of  the  primitive  variables  and  expresses  the  governing  equations 
in  conservation  form: 
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(pu)„  + (pv)  = 0, 

^ y 

(4) 

(pu^  + P)^  + (puv)  = 0, 

(5) 

(puv)  + (pv^  ♦ p)  = 0, 

A y 

(6) 

P = k . p"^. 

(7) 

Equations  (4)  to  (6)  all  have  the  same  structure  and  are  linear  in  the 
groups  of  variables.  In  the  present  formulation  advantage  is  taken  of  these 
features  to  reduce  substantially  the  number  of  cross-terms  that  appear. 
Equations  (4)  to  (6)  and  the  energy  equation. 
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have  been  found  to  be  suitable  for  solving  flow  problems  at  low  to  inter- 
mediate free-stream  Mach  numbers.  The  primitive  variable  formulation  gives 
the  required  final  solution,  the  pressure  on  the  body,  directly.  A dis- 
advantage of  the  primitive  variable  formulation  is  that  three  unknowns  per 
node  are  required. 

2.2  Equations  of  motion  and  finite  element  representation 

Because  of  the  technique  used  to  provide  starting  data  for  the  Galerkin 
formulation  (Section  3.3)  the  unsteady,  compressible,  inviscid  equations  will 
be  presented.  Following  Peyret  and  Viviand(ref .5)  the  conservation  form  of 
the  governing  equations  is 


aw  3f  dG 

at  ^ ax  ^ Sy 


0, 


(9) 


where  W,  F and  G are  three  component  vectors 


PU 

Pu^  + p 

puv 

PV 

. F = 

Puv 

. G = 

pv^  + p 

_p  _ 

-Pu 

_ PV 

W = 

Equations  (4)  to  (6)  are  equivalent  to 

ap  dG 


(10) 


= 0- 


(11) 


To  stabilise  the  integration  of  the  unsteady  equations  (Section  3.3)  the 
following  related  equations  will  also  be  considered, 


^ If  . 3G 
3t  * 3x 


a^w  aiw 

’[_3  x^  * 3y 


(12) 


The  coefficients  a,  in  equation  (12),  are  chosen  to  be  as  small  as  possible 
consistent  with  the  convergence  of  the  integration  of  the  unsteady  equations. 
Equations  (4)  to  (12)  are  non-dimensionalised  by  defining 


nd 


u/U^  %d  = 


■oo.  Pnd  = Pnd 


(P  - PcJ/PooUL-  (13) 


The  resultant  form  of  equations  (12)  remains  unaltered  if  a new  t and  a are 
defined.  From  now  on  the  subscript  nd  will  be  dropped.  Equation  (7) 
becomes 


1 + 7 . . p = p' 

and  equation  (8),  with  some  rearrangement,  becomes 


(14) 


1+7 


ft 


p = P 1 + 


7 - 1 


1 - (u'  + v^ ) 


(15) 


In  this  report  either  equation  (14)  or  (15)  has  been  used  to  link  p to  the 
other  variables. 
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The  far-field  boundary  conditions  are  applied  at  a finite  distance  from 
the  body  in  the  form 


u = u 


fs 


fs 


P = P 


fs’ 


(16) 


Two  choices  for  the  far-field  boundary  conditions  have  been  considered. 

Firstly  u^^  = 1.0,  = 0,  p^^  = 1.0.  Secondly  a Prandtl-Glauert 

transformation  has  been  applied  and  the  distorted  body  shape  has  been  used. 

For  the  flow  about  circular  cylinders  and  ellipses  the  far-field  velocity 
components  have  been  calculated  from  a complex  variable  solution  of  the 
distorted  body.  For  the  flow  about  aerofoils,  thin-aerofoil  theory  has  been 
applied  to  the  distorted  body  (see  Appendix  I).  Once  the  velocity  components 
are  available  equations  (14)  and(15)  have  been  used  to  give  the  density. 

Complex  variable  theory  and  thin  aerofoil  theory  have  also  been  used  to  give 
starting  data  throughout  the  flow- fie  Id. 

The  boundary  condition  at  the  body  has  required  that  the  flow  is  locally 
tangential  to  the  body  surface.  This  has  given  a relationship  between 
Uj  and  Vj  and  (/^)^  and  (^)j  at  the  nodes  on  the  body. 

As  the  first  step  of  the  finite  element  formulation  analytic  representations 
for  the  dependent  variables  are  introduced.  In  the  present  report  these  are 
introduced  for  the  groups  of  variables  W,  F and  G in  equation  (10).  E.g. 


pu  = 


Nj (x,y) 


(17) 


1 


• th 


node  and  indicates 


where  is  the  shape  function  appropriate  to  the  j 

the  nodal  value  of  the  different  variables.  It  is  believed  that  this  is 
the  first  finite  element  formulation  in  which  groups  of  variables,  rather  than 
single  variables,  have  been  given  an  analytic  representation.  This  technique 
has  been  used  Previously(ref.6,7)  in  the  application  of  an  Orthonormal  Method 
of  Integral  Relations  to  supersonic  boundary- layer  flow.  An  immediate 
advantage  of  finite  element  representations  like  equation  (17)  is  that,  since 
the  governing  equations  are  linear  in  the  groups  of  variables,  only  single 
summations  occur  after  application  of  the  finite  element  method.  ITius  the 
computation  of  the  equation  residuals  can  be  accomplished  more  efficiently. 

One  difficulty  associated  with  the  finite  element  method  applied  to  a 
system  of  equations,  like  (11)  and  (12),  is  to  ensure  that  the  orders  of  the 
shape  functions  that  appear  in  the  analytic  representations  like  equation  (17) 
are  chosen  consistently.  If  the  Galerkin  formulation  is  used  then  the 
corresponding  choice  of  the  order  of  the  weight  functions  presents  similar 
difficulties.  Failure  to  choose  the  orders  of  the  shape  functions 
consistently  produces  a less  efficient  solution  i.e.  more  nodal  unknowns  and 
a more  refined  grid  will  be  required  to  achieve  comparable  accuracy. 

Taylor  and  Hood(ref.8),  after  applying  a Galerkin  finite  element 
formulation  to  slow  viscous  flow,  concluded  that  for  consistency: 


(i)  the  maximum  order  of  error  associated  with  the  residual  of  each 
variable  must  be  equal 

(ii)  The  residuals  arising  from  each  equation  must  be  weighted  according 
to  the  maximum  error  occurring  in  each  equation. 
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It  is  well-known  that  for  slow  viscous  flow  the  momentum  equations  are 
dominated  by  the  pressure  gradient  and  viscous  dissipation  terms  whereas  the 
convective  or  inertia  terms  are  of  smaller  magnitude.  However  for  high  speed 
flows  the  momentum  equations  are  dominated  by  the  balance  between  the  inertia 
terms  and  the  pressure  gradient  terms;  the  viscous  terms  only  become 
significant  where  a large  velocity  gradient  occurs  e.g.  adjacent  to  a surface. 
Thus  it  would  seem  desirable  to  replace  (i)  with  the  condition  that 

(iii)  the  order  of  the  representation  of  the  flow  variables  should 

produce  errors  consistent  with  the  physical  processes  being  modelled. 

For  the  case  of  slow  viscous  flow  condition  (iii)  leads  to  a first  order 
representation  for  p and  a second  order  representation  for  u,  v. 

Equations  (4)  to  (6)  could  be  applied  to  incompressible,  inviscid  flow. 

In  this  case  p is  constant.  On  physical  grounds  it  would  be  expected  that 
all  groups  of  variables  in  each  equation  would  be  of  the  same  order  of 
magnitude.  Since  only  first  derivatives  appear  each  group  of  variables  in 
each  equation  will  require  the  same  order  of  representation  to  satisfy 
condition  (iii). 

Thus 


X = 


(x.y) 


j 


'18) 


where 


and 


X = (pu,  pv) 


Y = 


j 


(19) 


where  Y = (pu^ , puv,  pv^ , p) . L.  and  M.  are  shape  functions  of,  as  yet, 

^ ^ th 

undetermined  orders.  Suppose  u,  v are  represented  by  n order  shape 
functions  i.e. 


u 


n 


X 


then 


pu 


n 


X 


and 


p ~ Pu* 


2n 


X 


Thus  the  lowest  consistent  representation  would  be  for  L to  be  first  order 
and  M to  be  second  order,  i.e. 
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In  contrast  to  slow  viscous  flow,  inviscid  incompressible  flow  requires  a 
first  order  representation  for  u and  v and  a second  order  representation  for 
P- 

If  equations  (4)  to  (6)  are  applied  to  compressible,  inviscid  flow 
representations  (18)  and  (19)  are  also  applicable.  In  this  case 


P ~ 


P 
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therefore 
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7-1 


~ u 
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and  if 


u ~ X 
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Pu  ~ X 


* 1 
^ - 1 


) n 


and 


p ~ pu^ 


and 


O(M^) 

0(Lj) 


27 

7 + r 


\ 

(21) 


Since  only  integral  order  representations  are  possible  the  lowest,  consistent 
representation,  with  7 = 1.4,  would  have  as  a sixth  order  shape  function 

ard  Mj  as  a seventh  order  shape  function. 


If  a Galerkin  formulation  were 


employed  seventh  and  sixth  order  shape  functions  weighting  the  continuity 
and  momentum  equations  respectively  would  also  be  required.  Clearly  such  a 
high-order  representation  would  be  unwieldy  and  possibly,  because  of  the 
relatively  dense  stiffness  matrix,  be  inefficient. 

If  compressible,  inviscid  flow  is  to  be  represented  by  a moderately  low 
order  system,  equation  (21)  suggests  that  and  should  be  of  the  same 

order.  Consequently,  based  on  the  results  of  reference  9 all  groups  of 
variables  have  been  represented  by  quadratic  shape  functions  of  the  Serendipity 
family  in  rectangular  isoparametric  elements. 
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2.3  The  Galerkin  finite  element  formulation 
The  Galerkin  method  may  be  written  as 


N.  . R 
1 


(k) 


. dx  dy  = 0,  k = Ij  3; 


= 1,  n 


(22) 


where  is  the  shape  function  appropriate  to  the  i^^  node  and  R^*^^  is  the 

equation  residual  after  substitution  of  analytic  representations  like 
equation  (17).  For  the  present  problem  the  Galerkin  formulation  follows 
reference  10.  After  application  of  the  Galerkin  formulation  equation  (11) 
becomes 


dW. 


a.  . 
ij 


F.  + ; b. . . G. 

ij  1 


o / c..  .W.,i=l,n 


where 


(23) 


and 


s.  . 
ij 


a.  . 
11 


c.  . 
11 


1 1 dx  dy. 


dN. 

Ni  . . dx  dy. 

dN. 

Ni  . -5^  • dx  dy. 


dN.  dN. 

N.  . -T— ^ . dy  - -r- ^ • dx 
1 ox  ■'ay 


Ldx  ■ dx  ^ d 


dN.  dN. 

1 


dNj^ 

dy  - 


dx  dy. 


(24) 


. th 


node. 


Wj  etc.  are  the  values  of  W etc.,  given  by  equation  (12),  at  the  j 

The  coefficients  s. a. .,  b. . and  c. . are  evaluated  for  each  i^^  node  after 
11  11  11  11 

introducing  an  isoparametric  formulation (ref .9) . s^ ^ , a^^  and  b^j  can  be 

obtained  from  intermediate  coefficients  that  are  evaluated,  once  and  for 
all,  on  a dummy  element(ref .9) . a^^  and  b^^  are  the  same  expressions  as 

arise  in  the  treatment  of  inviscid,  incompressible  flow(ref.9). 
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2.4  Least-squares  finite  element  formulation 

The  least-squares  formulation  presented  here  follows  that  of  Zienkiewicz, 
Owen  and  Lee(ref.ll).  Substitution  of  the  analytic  representation, 
equation  (17),  into  the  governing  equations  (4)  to  (6),  produces  residuals 
of  the  following  form  in  each  element: 


• CD 


9N.  dN. 

j 


(25) 


• (2) 


dN.  _ _ dN. 

^ . (PU^  " ^ 


3 


(Puv)^, 


(26) 


• CD  . 


dN . dN  . 

^ . (PUV)  . . ^ • CPV^  ^ P). 


C27) 


The  least-squares  formulation  requires  that 


[ (tti  . + 02  . + aj  . R^^^^)  dx  dy  = minimum. 


C28) 


where  Oi  , 02  and  03  are  scalars  that  may  be  used  to  adjust  the  weight  of  the 
various  equations.  Differentiating  equation  (28)  with  respect  to  each  of  the 
unknown  nodal  values  in  turn  produces  the  following  result:- 

a.  . ■||‘‘’  . r(')  . a,  . If"’  . . a,  . dx.  dy  - 0 , i - 

(29) 


where  q^  = |(/fu)^,  (Pv)^,  Pj^i  • Substitution  of  equations  (25)  to  (27)  into 

equations  (29)  and  evaluation  of  the  integrals,  produces  the  following 
algebraic  equations: 


S.' 

1 


(m)  _ 


rf™^  . (pu) . + sf™^  . (pv)  . + tf™^  . (Pu^ ) . * . (Puv) . + 

ij  ^ D ij  ^ 3 ij  ^ 3 iJ  3 


(m)  -—2  > (m)  — 

y:  . . (pv  ) . ♦ z; . . p . 

D i]  ^3 


= 0,  m = 1,  3;  i = l,n. 


(30) 
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In  equation  (30)  m = 1 corresponds  to  in  equation  (29).  Thus 

three  equations  are  formed  at  each  node  if  there  are  three  unknown  nodal 
values.  In  equation  (30)  rf?^  etc.  are  algebraic  functions  of 
c. d. . , pU. , Wv.  and  p.  where 

ij  ij  i’  1 1 


a.  . 
ij 


c.  . 
ij 


d.  . 
ij 


//^ 

dN. 

• ^ 

3N. 

!!w 

■ ^ 

ff 

dN. 

//^ 

• 

ff 

dN. 

II W 

• ^ 

dx  dy. 


. dx  dy. 


. dx  dy , 


. dx  dy. 


y 


(31) 


Equations  (31)  may  be  compared  with  equations  (24)  which  arise  through 
application  of  the  Galerkin  formulation.  Once  isoparametric  elements  are 
introduced(ref . 9) , equations  (31)  are  considerably  more  complicated  and 
time-consuming  to  evaluate  than  equations  (24) . The  algebraic  expressions 

depend  on  whether  equation  (14)  or  (15)  is  used  to  obtain 


for  r. 


etc. 


p as  a function  of  the  other  variables, 
in  Appendix  II. 


The  detailed  expressions  are  given 


3.  ITERATIVE  SOLUTION  TECHNIQUES 

The  main  difficulty,  associated  with  solving  either  equations  (23)  that  arise 
from  the  Galerkin  formulation  or  equations  (30)  that  arise  from  the  least-squares 
formulation,  is  that  the  nodal  unknowns  /JUj , and  p^  occur  nonlinearly. 

Groups  of  terms,  like  puv^  , are  interpreted  as  (^^  . pvp/Wy 

Newton's  method  (Section  3.1)  was  used  initially  to  solve  the  steady  version 
of  equations  (23)  in  the  form 


a.  . 
ij 


F.  + 
J 


b.  . 
IJ 


. G. 
J 


= 0,  i = l,n. 


(32) 


However  this  produced  a singular  Jacobian.  This  problem  was  overcome  by 
differentiating  the  energy  equation  (15)  with  respect  to  x and  y and  using  these 
equations  instead  of  the  x-momentum  equation  on  the  y-axis  of  symmetry  and 
instead  of  both  momentum  equations  at  the  body  surface. 

In  practice  it  was  found  that  the  inversion  of  the  Jacobian  required  large 
amounts  of  main  storage  and  large  execution  time  even  after  introduction  of  the 
sparse  matrix  techniques  described  in  reference  9.  In  order  to  make  Newton's 
method  more  efficient  various  modifications  were  made  and  these  are  described 
in  Section  3.2.  As  the  number  of  nodal  unknowns  in  the  flow  field  was  increased 
difficulty  occurred  in  obtaining  starting  solutions  that  were  sufficiently  close 
to  the  converged  solution  to  permit  convergence  of  Newton's  method.  This  result 


J 


I 

i 
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is  consistent  with  the  findings  of  reference  12. 

To  alleviate  the  situation  the  governing  equations  in  unsteady  form  (23)  were 
integrated  as  though  they  were  simultaneous  ordinary  differential  equations  in 
time.  This  was  a fairly  slow  process  and  was  terminated  as  soon  as  the  iterative 
solution  was  close  enough  to  the  final  solution  to  permit  Newton's  method  to 
converge.  A description  is  given  in  Section  3.3.  Because  the  least-squares 
formulation  leads  to  a positive-definite  stiffness  matrix  it  is  possible  to  use 
an  SOR-type  iterative  scheme  to  solve  equations  (30).  This  is  described  in 
Section  3.4. 

3.1  Newton's  method 

This  will  be  described  in  relation  to  the  solution  of  equation  (32)  which 
can  be  written  as 


Ri(q)  = 0,  i = l.n.  (33) 

where  q is  a vector  of  all  the  nodal  unknowns  w.,  equation  (10),  and  R^ 

represents_the  residual  of  the  i'^"  equation.  If  an  arbitrary  starting 
solution,  q^,  is  substituted  into  equation  (32)  non-zero  equation  residuals 

result.  In  the  neighbourhood  of  q^  the  solution  can  be  obtained  by 

application  of  a Taylor  series  expansion. 

= ^o  ^ If 

If  the  Taylor  series  is  truncated  as  shown  and  if  it  is  assumed  that  qi  is 
the  exact  solution  then  Rj  = 0 and 

q.  = % - J"'  (%)  . 

where  J is  the  Jacobian,  3R/3q.  If  equation  (35)  is  interpreted  as  one  step 
of  an  iteration  Newton's  method  is  obtained: 


‘'(^+1  = % - 

The  effectiveness  of  Newton's  method  depends  on  the  accuracy  of  the 
assumptions  underlying  equations  (34)  and  (35).  Clearly  if  the  starting 
solution  is  close  to  the  converged  solution  the  assumptions  are  reasonable 
and  Newton's  method  is  convergent;  close  to  the  converged  solution  the 
method  has  the  property  of  quadratic  convergence(ref . 13) . A discussion  of 
Newton's  method  and  some  iterative  techniques  based  on  Newton's  method  are 
given  in  reference  14. 

3.2  Modified  Newton's  method 

In  practical  applications  most  of  the  execution  time  is  spent  factoring 

to  form  J~'  even  if  sparse  matrix  techniques  and  the  storage  of  J~'  in 

factored  form  are  used.  Typically,  in  early  applications  of  the  Galerkin 
formulation  to  the  present  problem,  the  factorisation  of  the  Jacobian 
accounted  for  80  to  90%  of  the  execution  time.  Any  modification  that 
permits  more  utilisation  of  each  evaluation  of  J”'  is  clearly  desirable. 

An  obvious  modification  is  to  compute  p steps  with  the  same  J^  i.e. 
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% -f  1 ~ % 


0 < V < p - 1. 


(37) 


Another  improvement  can  be  obtained  if  is  treated  as  a vector  6^^ 

which  is  small  compared  with  . Then  6^  may  be  interpreted  as  a search 

direction.  A series  of  residuals  Rj,  are  evaluated  for  various  correspond- 

. — X , 

ing  to  given  by 


-k 

% + l 


% - \ • 'r- 


(38) 


Xj^  is  a scalar  and  typically  takes  values  between  0 and  2.  If  ^ 
corresponds  to  the  solution  q^  i.e.  X^  = 0,  then  two  additional  solutions, 
RJ,  and  R^  are  obtained  corresponding  to  Xi  , Xj  . For  each  solution  the 

global  sum  of  all  the  residuals,  , is  computed  from 


(39) 


A quadratic  dependence  of  on  X^^  is  extracted  and  the  value  X^^^,  which 

minimises  F^^  is  used  to  obtain  The  technique  of  finding  a minimum  in 

a particular  search  direction  has  been  used  previously  for  minimising  a sum 
of  squares (ref . 15) . 

Another  modification  to  equation  (36)  is  possible  by  attempting  to 
approximate  J ' . Assuming  that  J is  close  to  J , let 


= J + 7. 
o 


If  a first  approximation  to  is  confuted  as 


= J*'  R , 
o u • 


then  a better  approximation  is 

8 


- 


(40) 


(41) 


(42) 


where  is  an  approximation  to  • From  equation  (40) 


= j„(i  - j;'  T) 


and 


= (I  + J 


7)'' 


(^3) 


. J’‘ 
o 


r- 
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Equation  (43)  can  be  approximated  by 


J.'k  = - '^o  • '^o  • 


Substitution  into  equation  (42)  gives 


(I  - j;-  . 7)  . j;'  . R, 


(I  - j;'  . 7)  . 6/. 


Substitution  for  7 and  simplification  gives 

= V - - V • V 

= 2^*  - j;‘  . . 6/.  (44) 

The  modified  Newton's  method  used  in  the  present  study  has  consisted  of  the 
following  steps: 

(i)  compute  and  J^' 

(ii)  use  equation  (42)  to  compute  6^ 

(iii)  find  a minimum  F in  the  6 direction 
_ o o 

(iv)  compute  q^^ . Jj, 

(v)  use  equation  (41)  to  compute  6^* 

(vi)  use  equation  (44)  to  compute  6^ 

(vii)  compute  a minimum  in  the  direction 
(viii)  if  Fj^  not  sufficiently  small  go  to  step  (iv) . 

3.3  Starting  data  for  Newton's  method 

Equation  (23)  can  be  rewritten  as 


V Ft  , 

‘ii  • ^ = 


- •'i'  i 


where 


R.  = 7 a...F.+7  b...G.-a7  c...W.. 

1 Z_j  ij  J Z_y  ij  J Z_/  ij  J 


In  matrix  form  equation  (45)  can  be  inverted  to  give 


W = - S . R. 
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The  coefficients  of  the  matrices  S only  depend  on  the  shape  functions  and 
consequently  can  be  inverted  once  and  for  all.  To  avoid  excessive  fill-in 
and  to  speed  up  the  matrix  multiplications  on  the  right  hand  side  of  equation 
(46),  matrices  like  S'*  are  stored  in  sparse  factored  form. 

If  equations  (46)  are  integrated  for  large  time,  the  residuals  approach 

zero  and  consequently  the  time  derivatives  also  approach  zero.  The 

different  R^^  approach  zero  at  significantly  different  rates  and  this  causes 

difficulties  for  the  efficient  use  of  the  integration  routine.  Gear(ref.l6) 
has  produced  a predictor/corrector  algorithm  to  suit  this  situation.  In 
addition  Gear's  algorithm  automatically  adjusts  the  step-size  and  the  order 
of  the  predictor/corrector  formula  to  minimise  the  integration  error. 

The  term  on  the  right  hand  side  of  equation  (23)  has  been  introduced  to 
stabilise  the  numerical  integration  of  equations  (46).  Although  the 
additional  terms  have  some  physical  basis (ref. 5)  for  the  momentum  equations 
there  is  no  such  basis  for  the  continuity  equations.  Therefore  it  is 
desirable  that  a should  be  as  small  as  possible  consistent  with  convergence. 

_ By  expanding  equation  (46)  as  a Taylor  series  about  a starting  solution 
q^  it  is  possible  to  carry  out  an  eigenvalue  analysis  to  determine  the 

critical  value  of  a,  below  which  the  integration  of  equation  (46)  is  unstable. 
The  critical  value  was  found  to  be  approximately  a = 0.1,  and  this  was 
confirmed  in  actual  computations. 

3.4  One  step  SOR-Newton's  method 


Application  of  the  least-squares  finite  element  formulation  (Section  2.4) 
produces  a positive  definite  stiffness  matrix.  Consequently  an  SOR-type 
solution  technique  is  possible.  This  may  be  obtained  by  considering  a small 

change  in  (equation  (30))  due  to  a small  change  in  q^,  where  q^  is  either 

/nr,  /nr.,  p.. 


s(m) 


r . 

(q  ) - 

1 '■^o  oq^ 


Aq.  = 


0, 


(47) 


or  in  a more  general  form 


v + 1 
^i 


‘li 


X 


as.^'"^(q  ) 

L aq.  J 


(48) 


Examination  of  equation  (30)  and  Appendix  II  indicates  that  both  terms  like 

(pu^)j  and  etc.  are  functions  of  q^ , thus  SS^^V^q^  is  algebraicly 

complicated.  In  contrast  to  the  full  Newton's  method  3s^/3q^  is  a scalar 

and  trivial  to  invert,  thus  no  excessive  demand  is  made  on  main  storage  or 
computation  time.  However  the  rate  of  convergence  is  slower  than  Newton's 
method.  A scalar,  X,  has  been  used  to  increase  the  rate  of  convergence. 

If  X was  greater  than  1.7  the  iterative  process  diverged. 


4.  FLOW  ABOUT  CIRCULAR  AND  ELLIPTIC  CYLINDERS 

Numerical  solutions  of  the  flow-field  about  circular  cylinders  have  been 
obtained  using  a Galerkin  finite  element  formulation  and  a least-squares  finite 
element  formulation.  A schematic  representation  of  the  grid  used  to  obtain 
solutions  about  both  circular  and  elliptic  cylinders  is  shown  in  figure  1.  The 
grid  system  is  essentially  polar.  An  isoparametric  formulation (ref .9)  has  been 
used  to  relate  this  to  a cartesian  coordinate  system. 
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Two  stagnation  points  exist  for  the  flow  about  bluff  bodies.  At  the  stag- 
nation points  the  velocity  components,  u v are  zero.  In  the  present  formulation 
the  pressure  and  density  at  the  stagnation  points  have  been  considered  part  of 
the  boundary  conditions  by  assuming  a streamline  attaches  both  stagnation  points 
to  the  freestream.  Consequently  equations  (14)  and  (15)  can  be  used  to  give 
the  following  stagnation  point  values  for  p and  p, 


•i 


and 


(49) 


sp 


(50) 


The  far-field  boundary  conditions  are  obtained  by  scaling  the  x-coordinate  by 
the  Prandtl-Glauert  factor  and  then  computing  the  far-field  velocity  components 
using  complex  variable  theory  (Appendix  I). 

4.1  Circular  cylinder 

Results  are  presented  in  figure  2 for  the  variation  of  the  surface  Mach 
number  with  angular  location.  Results  for  free-stream  Mach  numbers  0.1  to 
0.3  have  been  obtained  using  the  Galerkin  finite  element  formulation.  These 
results  have  been  obtained  with  91  elements  and  465  nodal  unknowns  spanning 
the  flow- field.  It  was  found  that  for  free-stream  Mach  numbers  greater  than 
0.32  Newton's  method  would  not  converge  even  when  starting  from  a converged 
solution  at  a free-stream  Mach  number  only  marginally  smaller.  The  reason 
for  this  is  not  known. 

Solutions  obtained  using  the  least-squares  formulation  at  freestream 
Mach  numbers  of  0.2  and  0.3  are  also  shown  in  figure  2.  These  results  were 

obtained  with  91  elements  and  829  nodal  unknowns  spanning  the  flow-field. 

It  is  apparent  that,  as  the  free-stream  Mach  number  is  increased,  the 
solutions  obtained  with  the  Galerkin  formulation  indicate  a significantly 
greater  acceleration  of  the  flow  than  those  obtained  with  the  least-squares 
formulation. 

Results  for  the  flow  about  a circular  cylinder  at  a free-stream  Mach 
number  of  0.4  are  presented  in  figure  3.  The  finite  element  solution  has 
been  obtained  with  a least-squares  formulation.  The  results  presented  for 
the  surface  Mach  number  have  been  obtained  with  91  elements  and  829  nodal 
unknowns  spanning  the  flow-field. 

The  finite  element  solution  is  compared  with  solutions  obtained  by  the 
Method  of  Lines (ref . 17)  and  series  solutions  presented  by  Greenspan  and 
Jain(ref . 18) ; the  series  solutions  are  due  to  Imai(ref.l9)  and  Lush  and 
Cherry(ref . 20) . The  solution  by  the  Method  of  Lines  agrees  closely  with 
those  of  references  19  and  20  and,  consequently,  has  not  been  plotted. 

The  finite  element  solution  is  in  reasonable  agreement  with  the  other 
solutions  although  it  underpredicts  the  other  solutions  close  to  the  cylinder 
shoulder  point.  Since  the  solutions  given  by  references  19  and  20  have  been 
obtained  by  a truncated  series  representation  the  difference  between  those 
solutions  and  the  finite  element  solution  is  not  considered  significant. 

4.2  Elliptic  cylinder 


A least-squares  finite  element  solution  for  the  flow  about  a 2:1  elliptic 
cylinder  at  a free-stream  Mach  number  of  0.5  is  shown  in  figure  4.  These 
results  were  obtained  with  91  elements  and  829  nodal  unknowns  spanning  the 
flow-field.  Solutions  were  also  obtained  using  the  Galerkin  formulation  at 
free-stream  Mach  numbers  up  to  0.3  but  convergence  did  not  occur  for  major 
axis/minor  axis  ratios  greater  than  1.8. 
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Numerical  calculations  using  the  Method  of  Integral  Relations(ref . 21) 
indicate  that  the  flow  is  just  supersonic  at  the  shoulder  of  the  elliptic 
cylinder  for  a free-stream  Mach  number  of  0.5.  In  contrast  the  present 

solutions  indicate  that  the  flow  is  everywhere  subcritical.  The  present 
results  are  compared  with  solutions  obtained  by  the  Method  of  Lines (ref . 17) 
which  are  shown  in  figure  4.  Agreement  is  quite  good  although,  as  with  the 
results  in  figure  3,  the  finite  element  solution  slightly  underpredicts  the 
surface  velocity  on  the  'plateau' . 


5.  FLOW  ABOUT  AEROFOILS 

Solutions  for  the  flow- field  about  two  representative  aerofoils  have  been 
obtained  using  the  least-squares  formulation.  The  nature  of  the  grids  used  has 
been  dictated  partially  by  the  shape  of  the  aerofoil  in  question.  A schematic  of 

the  grid  used  for  the  circular-arc  aerofoil  is  shown  in  figure  5.  A deliberate 

attempt  has  been  made  to  keep  elements  far  removed  from  the  body  surface 
rectangular;  all  elements  outside  the  region  ABCD  are  rectangular.  This  reduces 
the  number  of  cross-terms  and  hence  the  computation  time  required  to  manipulate 
those  terms . 

The  NACA-0012  aerofoil  has  a forward  stagnation  point  and  if  the  type  of  grid 
shown  in  figure  5 were  used  gross  distortion  of  the  elements  adjacent  to  the 
forward  stagnation  point  would  occur.  Although  the  isoparametric  formulation  is 
successful  in  mapping  rectangular  elements  onto  a non-rectangular  region,  the 
quality  of  the  solution  obtained  is  degraded  if  the  elements  become  distorted. 

A schematic  of  the  grid  used  to  obtain  flow-field  solutions  about  the  NACA-0012 
aerofoil  is  shown  in  figure  6.  All  elements  inside  the  area  ABCDEF  are 
rectangular. 

A comparison  of  the  grids  shown  in  figures  1,  5 and  6 indicates  that  basically 
two  types  of  grid  have  been  used  depending  on  the  local  body  geometry.  Where 
the  body  shape  can  be  locally  approximated  by  the  arc  of  a circle  a polar  grid  has 
been  used.  Where  the  body  shape  can  be  locally  approximated  by  a straight  line 
a cartesian  grid  has  been  used. 

The  local  values  p and  p at  the  forward  stagnation  point  of  the  NACA-0012 
aerofoil  are  given  by  equations  (49)  and  (50).  The  far- field  boundary  conditions 
for  both  aerofoils  have  been  obtained  by  applying  a Prandtl-Glauert  transformation 
to  the  x-coordinate  and  applying  thin-aerofoil  theory  to  the  distorted  body. 

This  is  described  in  Appendix  I. 

5.1  6%  circular-arc  aerofoil 

Solutions  to  the  flow  about  a 6%  circular-arc  aerofoil  have  been  obtained 
for  various  free-stream  Mach  numbers.  All  the  results  presented  in  figures  7 
to  10  have  been  obtained  with  102  elements  and  868  nodal  unknowns  spanning 
the  flow-field. 

The  surface  pressure  distribution  for  a free-stream  Mach  number  of  0.71  is 
shown  in  figure  7.  Also  shown  in  figure  7 are  experimental  results  due  to 
Knechtel (ref . 22) . Knechtel  obtained  pressure  distributions  for  a completely 
smooth  aerofoil  for  which  the  boundary  layer  would  have  been  laminar  at 
least  to  the  50%  chord  point.  By  introducing  roughness  just  aft  of  the 
leading  edge  Knechtel  also  obtained  results  for  which  the  boundary  layer  was 
turbulent  throughout.  An  examination  of  figure  7 indicates  that  the  finite 
element  solution  is  in  good  agreement  with  the  experimental  solution  obtained 
with  a smooth  aerofoil.  The  shape  of  the  pressure  distribution  is  slightly 
different  particularly  close  to  the  leading  and  trailing  edges.  This  may  be 
due  to  the  courseness  of  the  grid  at  the  body  surface  used  to  obtain  the 
computational  solution.  A consideration  of  the  experimental  accuracy 
suggests  that  in  an  area  of  rapidly  changing  pressure  such  differences 
indicated  may  not  be  significant. 

The  surface  pressure  distribution  for  a free-stream  Mach  number  of  0.82  is 
shown  in  figure  8.  Also  included  in  figure  8 are  experimental  results  due  to 
Knechtel (ref . 22) . The  experimental  results  were  obtained  with  a smooth 
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aerofoil.  As  in  figure  7 the  agreement  shown  in  figure  8 is  good.  For  the 
results  shown  in  figures  7 and  8 the  flow-field  is  everywhere  subsonic. 

The  results  shown  in  figures  9 and  10  have  been  obtained  at  a free-stream 
Mach  number  of  0.88  for  which  the  surface  flow  is  just  sonic  at  the  thickest 
part  of  the  aerofoil.  The  pressure  distribution  for  this  case  is  shown  in 

figure  9.  Experimental  results  at  = 0.88  were  not  obtained  by  Knechtel 

(ref. 22),  although  an  interpolation  of  the  results  that  are  available 
suggests  that  the  negative  pressure  coefficient  at  the  50%  chord  point  is 
larger  for  the  experimental  results  than  it  is  for  the  computational  results 
shown  in  figure  9.  It  is  likely  that  the  displacement  thickness  effect  has 
caused  some  supersonic  flow  and  this  has  increased  the  suction  peak.  The 
corresponding  surface  Mach  number  distribution  is  shown  in  figure  10  and  the 
local  Mach  number  at  the  50%  chord  point  can  be  seen  to  be  sonic. 

5.2  NACA-0012  aerofoil 

Computational  solutions  have  been  obtained  for  the  flow  about  a NACA-0012 
aerofoil  at  free-stream  Mach  numbers  of  0.4  and  0.72.  The  finite  element 
solution  for  the  free-stream  Mach  number  of  0.4  has  been  obtained  with  a least- 
squares  formulation  and  88  elements  and  767  nodal  unknowns  spanning  the  flow- 
field.  This  solution  is  shown  in  figure  11.  Also  shown  in  figure  11  are 
experimental  results  due  to  Amick (ref . 23)  and  some  computational  results, 
using  a finite  difference  method,  due  to  Emmons (ref . 24) . 

The  finite  element  solutions  have  simulated  an  aerofoil  in  an  unconstrained 
free-stream.  In  contrast  the  experimental  results  of  Amick  have  not  been 
corrected  for  the  influence  of  the  wind-tunnel  walls  and  the  computational 
results  of  Emmons  have  deliberately  allowed  for  the  presence  of  the  wind- 
tunnel  walls.  The  discussion  in  reference  24  suggests  the  experimental 
results  of  Amick  would  produce  a maximum  negative  pressure  coefficient 
approximately  4 to  5%  less  if  the  aerofoil  were  in  an  unconstrained  free-stream. 

It  is  apparent  that  the  occurrence  of  the  stagnation  point  in  the  flow- 
field  and  the  consequent  acceleration  of  the  flow  past  the  nose  of  the  aero- 
foil causes  a very  rapid  change  in  the  values  of  the  dependent  variables  in 
the  nose  region.  To  obtain  accurate  solutions  this  situation  requires  small 
elements  in  that  region  and,  if  the  computation  time  is  to  be  kept  within 
reasonable  bounds  requires  large  elements  elsewhere(figure  6).  The  automatic 
mesh  generation  scheme  is  the  same  as  that  used  and  described  in  reference  9 
and  the  above  requirements  have  only  been  partially  achieved.  Nevertheless 
the  agreement  with  the  experimental  results  shown  in  figure  11  is  good 
particularly  away  from  the  nose  region. 

A solution  obtained  for  the  flow  about  a NACA-0012  aerofoil  at  a free- 
stream  Mach  number  of  0.72  is  presented  in  figure  12.  These  results  were 
obtained  with  the  least-squares  formulation  and  96  elements  and  837  nodal 
unknowns  spanning  the  flow-field.  Computational  result  . for  this  problem 
have  been  obtained  by  Lock (ref. 25)  using  the  method  of  Sells(ref.26) . Lock 
considers  that  Sells'  method  is  capable  of  giving  results  that  are  accurate 
to  1%  of  the  maximum  perturbation  velocity,  A feature  of  Sells'  method  is 
the  mapping  of  the  flow-field  onto  the  interior  of  a unit  circle.  This 
avoids  the  problem  of  applying  the  far-field  boundary  condition  at  a finite 
distance  from  the  body.  The  solution  obtained  by  Sells'  method  is  shown 
in  figure  12.  An  examination  of  figure  12  indicates  that  the  finite  element 
solution  for  the  pressure  coefficient  underpredicts  Sells'  solution 
particularly  in  the  nose  region.  This  may  be  due  to  not  applying  the  far- 
field  boundary  conditions  sufficiently  far  from  the  body.  A more  likely 
cause  is  the  courseness  of  the  grid  in  a region  where  large  gradients  are 
occurring . 


kL.. 


17 


WRE-TR-1858(W) 


6.  COMPARISON  OF  THE  GALERKIN  AND  LEAST-SQUARES  FORMULATIONS 

The  Galerkin  finite  element  formulation  has  been  used  frequently  to  find  the 
numerical  solution  to  flow  problems.  However,  most  previous  applications  of  the 
Galerkin  formulation  have  been  to  the  full  Navier-Stokes  equations  in  incompress- 
ible form.  In  such  a situation  the  viscous  terms  in  the  governing  equations 
play  a stabilising  role  in  both  the  physical  and  numerical  sense.  Viscosity 
plays  no  role  in  the  present  problem  unless  artificially  introduced  (Section  3.3). 
Most  previous  finite  element  solutions  of  the  present  problem  have  been  based  on 
a variational  formulation  which  leads  to  a stiffness  matrix  that  is  positive 
definite . 

The  application  of  a Galerkin  finite  element  formulation  to  the  present  problem 
has  led  to  algebraic  governing  equations  that  are  algebraicly  simple,  economical 
to  create  and  relatively  sparse.  In  fact  the  derivation  of  the  coefficients 
^ij’  ^i j ’ (24),  in  the  algebraic  equations  is  precisely  the  same  as  that 

for  the  application  of  the  Galerkin  formulation  to  incompressible,  inviscid  flow 
(ref .9) . 

A problem  with  the  Galerkin  formulation  arises  in  relation  to  finding  a 
suitable  iterative  scheme  for  solving  the  nonlinear  algebraic  governing  equations. 
In  the  application  of  the  Galerkin  finite  element  method  to  slow  viscous  flow 
Newton's  method  has  been  used  with  considerable  success.  Since  Newton's  method 
is  completely  general  it  should  be  effective  in  the  present  situation.  The  first 
difficulty  that  occurred  in  implementing  Newton's  method  was  that  the  Jacobian 
was  singular  due  to  the  symmetry  about  the  y-axis.  To  overcome  this  problem  it 
was  necessary  to  introduce  special  equations  on  the  y-axis  and  on  the  body.  A 

second  difficulty  with  Newton's  method,  and  it  may  be  related  to  the  first 
difficulty,  was  that  it  failed  to  converge  for  large  numbers  of  nodal  unknowns 
even  when  starting  from  a solution  close  to  a previously  converged  solution. 

In  contrast  the  least-squares  finite  element  formulation  produced  algebraic 
equations  that  were  more  complex,  required  more  computation  time  to  form  and  were 
less  sparse  than  those  produced  by  the  Galerkin  formulation.  However,  because 
of  the  positive  definite  nature  of  the  stiffness  matrix  the  iterative  technique 
for  solving  the  algebraic  governing  equations  produced  no  convergence  problems. 

At  intermediate  stages  of  the  iterative  process  the  solution  was  quite  smooth 
which  was  in  marked  contrast  to  that  produced  by  the  Galerkin  formulation  at 
intermediate  stages. 

It  seems  reasonable  to  conclude  that  since  the  current  problem  is  dominated  by 
the  non-linear  convective  terms  the  choice  of  the  iterative  scheme  for  solving 
the  non-linear  algebraic  governing  equations  becomes  crucial.  Thus  any  finite 
element  formulation  that  leads  to  a positive  definite  stiffness  matrix  is  likely 
to  have  a considerable  advantage  over  any  formulation  that  does  not. 


7.  CONCLUSIONS 

The  main  conclusion  of  this  report  is  that  the  least-squares  finite  element 
formulation  has  been  very  effective  in  obtaining  solutions  to  compressible, 
inviscid  flow  and  the  Galerkin  formulation  has  not.  Additional  conclusions  are 
that  the  use  of  a primitive  variable  formulation  results  in  very  accurate 
solutions  and  that  the  combination  of  using  the  conservation  form  of  the  governing 
equations  and  representing  groups  of  variables  rather  than  single  variables  leads 
to  a sparser  stiffness  matrix  and  hence  a more  economical  solution.  In  contrast 
to  the  treatment  of  slow,  viscous  flow,  the  most  efficient  finite  element 
solution  of  the  present  problem  is  obtained  if  each  group  of  variables  in  the 
governing  equations  is  represented  by  shape  functions  of  the  same  order. 
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NOTATION 

vectors  of  groups  of  terms  appearing  in  the 
equations  of  motion 

unit  matrix 

Jacobian 

shape  function 

Mach  number,  sliape  function 

Mach  number  at  body  surface 

shape  function 

equation  residual 

free-stream  velocity 

vector  of  nodal  unknowns 

vector  of  groups  of  terms  appearing  in  the 
continuity  equation 

vector  of  groups  of  terms  appearing  in  the 
momentum  equations 

sound  speed 

algebraic  coefficients  in  the  governing 
equations  (Galerkin  formulation) 

order  of  shape  function,  number  of  nodes 

pressure 

vector  of  nodal  unknowns 

tangential  velocity  at  body  surface 

algebraic  groups  in  t .e  governing  equations 
(least-squares  formulation) 

coordinate  along  the  body  surface 

time 

velocity  component  in  the  x direction 
velocity  component  in  the  y direction 
cartesian  coordinates 


a 


parameter  controlling  stabilising  terms  in  the 
equations  of  motion;  parameter  controlling 
relative  influence  of  the  residuals 
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specific  heat  ratio;  correction  to  the 
Jacobian 

correction  vector 

scalar  in  Newton's  method 

density 

velocity  potential 

free-strcam 
nondimensional 
stagnation  point 
iteration  step 

conditions  in  the  free-stream 
nodal  value 

differentiation  with  respect  to  time 
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APPENDIX  I 


FAR-FIELD  BOUNDARY  CONDITIONS 

For  most  of  the  results  presented  in  this  report  the  far-field  boundary 
conditions  have  been  determined  by  assuming  that  the  body  introduces  a disturbance 
into  the  flow.  Since  the  disturbance  in  the  far-field  is  small  it  is  appropriate 
to  introduce  a disturbance  potential,  which  perturbs  a uniform  free-stream, 
i.e. , 


•h  = U^  . X + ip- 

For  subcritical  flow  ^ is  governed  by 

1 


(I.l) 


„ 

(1  - * 3y^ 


CI-2) 


A new  coordinate  x = (1  - . x is  introduced  and  the  governing  equation  (1.2) 

reduces  to  Laplace's  equation 


■ »• 


(1.3) 


For  the  flow  about  circular  and  elliptic  cylinders  complex  variable  solutions 
(which  satisfy  equation  (1.3))  are  introduced  and  the  far-field  velocity  compo- 
nents calculated. 

For  the  flow  about  aerofoils  a continuous  source  distribution  is  introduced 
along  the  chord  line  of  the  aerofoil.  The  source  strength  is  proportional  to 
the  local  distorted  body  slope,  dy/dx”.  The  source  distribution  satisfies 
equation  (1.3).  The  velocities  induced  by  the  source  distribution  are 
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s is  the  non-dimensional  distance  along  the  aerofoil  chord.  The  integration  in 
equation  (1.4)  has  been  performed  numerically  by  dividing  up  the  aerofoil  chord 
into  the  same  finite  elements  as  are  used  in  the  body  of  the  report,  and  summing 
the  contributions  from  each  element.  Gauss  quadrature  and  an  isoparametric 
formulation  are  utilised  to  perform  the  integration  over  each  element. 

Once  the  far-field  velocity  components  have  been  obtained  the  corresponding 
density  and  pressure  are  given  by 
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The  al)Ove  linearised  solutions  have  also  been  used  to  provide  starting  data 
throughout  the  domain. 
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APPENDIX  11 

DERIVATION  OF  THE  AI.GEBRAIC  EQUATIONS  OBTAINED 
FROM  THE  I, EAST-SQUARES  FORMULATION 

After  application  of  the  least-squares  finite  element  formulation  to 
equations  (4)  to  (6)  the  following  algebraic  equations  are  obtained: 

" ^ij^  • ^ ^i?  ■ * ^ij 


(pvM,  + . p,l  = 0 


J ij  " 1 


m = 1,3;  i = l,n. 


(II. 1) 


In  equations  (II.l)  the  actual  expressions  for  r^j  etc.  depend  on  which  relation- 
ship is  used  to  link  p to  the  other  variables.  If  p is  related  to  the  density 
by  the  isentropic  equation 


1 + 7 . P = 


then  the  following  expressions  are  obtained 
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If  p is  related  to  the  other  variables  through  the  energy  equation 

= P.|l  + ^4-^  • C I 1 - ♦ vMll 


the  following  expressions  for  r.^l*^  etc.  are  obtained 
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Figure  1. 


Flow-field  geometry  and  schematic  of  grid  for  the 
circular  cylinder 
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Figure  2.  Surface  Mach  number  variation  for  flow  about  a 
circular  cylinder 
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Figure  3.  Surface  velocity  variatioTi  for  the  flow  about  a circular 
cylinder  at  = 0.40 
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Figure  4.  Surface  velocity  variation  for  the  flow  about  a 2:1 
elliptic  cylinder  at  = 0.50 
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AEROFOIL  E 0 

Figure  6.  Schematic  of  the  finite  element  grid  used  for  the 
NACA-0012  aerofoil 
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